Eigenfunctions of the Koopman operator
Adapted from: Example 2.6 of [MAI20]
[MAI20] Mauroy, Alexandre, Aivar Sootla, and Igor Mezić. Koopman framework for global stability analysis. The Koopman Operator in Systems and Control: Concepts, Methodologies, and Applications (2020): 35-58.
using DynamicPolynomials
using LinearAlgebra
@polyvar x[1:2]
a = 1
I = 0.05
ε = 0.08
γ = 1
F0 = [-x[2] - x[1] * (x[1] - 1) * (x[1] - a) + I, ε * (x[1] - γ * x[2])]2-element Vector{DynamicPolynomials.Polynomial{Float64}}:
0.05 - x₂ - x₁ + 2.0x₁² - x₁³
-0.08x₂ + 0.08x₁We move the equilibrium, approximately (0.0256, 0.0256), to the origin. Since γ = 1, both coordinates are equal. Newton's method refines the rounded value so that the translated vector field vanishes at the origin.
x_eq = let z = 0.0256
for _ in 1:5
z -= (-z - z * (z - 1)^2 + I) / (-2 + 4z - 3z^2)
end
z
end
F = [f(x => x .+ x_eq) for f in F0]2-element Vector{DynamicPolynomials.Polynomial{Float64}}:
-3.7100043089738355e-18 - x₂ - 0.8993758546202244x₁ + 1.9230516279758776x₁² - x₁³
-0.08x₂ + 0.08x₁We compute the Jacobian at the equilibrium.
J = [j(x => zeros(2)) for j in differentiate(F, x)]2×2 Matrix{Float64}:
-0.899376 -1.0
0.08 -0.08Its eigenvalues have negative real part. The linear term w ⋅ x of a Koopman eigenfunction must satisfy J' * w = λ * w, so we need a left eigenvector of J. A right eigenvector leaves a nonzero linear residual which the higher-degree terms cannot cancel. We therefore use a left eigenvector instead of the right eigenvector stated in Example 2.6, consistently with Eq. (2.7) on p. 42 of the book.
E = eigen(J')
λ = E.values[end]
w = E.vectors[:, end]2-element Vector{Float64}:
-0.11258237486231153
-0.9936423948636461Eigenfunctions are defined up to a scalar factor. We choose a unit gradient with positive second component to fix the sign of the contour labels. Figure 2.3 appears to use the original coordinates, with the equilibrium at (x_eq, x_eq), and a different normalization. Our axes are shifted to the equilibrium and our gradient has unit norm, so the contour labels need not coincide with those in the book.
if w[2] < 0
w = -w
end2-element Vector{Float64}:
0.11258237486231153
0.9936423948636461We approximate the eigenfunction on the disk of radius r. Working in coordinates y = x / r avoids the very small high-degree monomials on this disk. The scaled vector field is G(y) = F(r * y) / r.
using SumOfSquares
r = 0.3
@polyvar y[1:2]
G = [f(x => r .* y) / r for f in F]
X = @set y[1]^2 + y[2]^2 ≤ 1Basic semialgebraic Set defined by no equality
1 inequality
1 - y[2]^2 - y[1]^2 ≥ 0
We define the program for the FitzHugh-Nagumo problem below: N is the degree of ϕN as defined in [MAI20, p. 50] and M is the degree of the multipliers for the constraints of X. As maxdegree corresponds to the degree of the multiplier multiplied by 1 - y[1]^2 - y[2]^2, we set it to M + 2.
function fitzhugh_nagumo(solver, N, M)
model = SOSModel(solver)
@variable(model, η)
@objective(model, Min, η)
@variable(model, ϕN, Poly(monomials(y, 2:N)))
ϕ = r * (w ⋅ y) + ϕN
residual = G ⋅ differentiate(ϕ, y) - λ * ϕ
@constraint(model, -η ≤ residual, domain = X, maxdegree = M + 2)
@constraint(model, residual ≤ η, domain = X, maxdegree = M + 2)
optimize!(model)
return ϕ(y => x ./ r), model
endfitzhugh_nagumo (generic function with 1 method)In [MAI20, p. 50], the result is obtained with N = 10 and M = 20.
import CSDP
solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)
ϕ, model = fitzhugh_nagumo(solver, 10, 20)
solution_summary(model)solution_summary(; result = 1, verbose = false)
├ solver_name : CSDP
├ Termination
│ ├ termination_status : ALMOST_OPTIMAL
│ ├ result_count : 1
│ └ raw_status : Problem solved to near optimality.
├ Solution (result = 1)
│ ├ primal_status : NEARLY_FEASIBLE_POINT
│ ├ dual_status : NEARLY_FEASIBLE_POINT
│ ├ objective_value : 7.98439e-06
│ └ dual_objective_value : 7.97639e-06
└ Work counters
└ solve_time (sec) : 4.96063e-01The optimal polynomial is obtained as follows:
ϕ_opt = value(ϕ)
residual_opt = F ⋅ differentiate(ϕ_opt, x) - λ * ϕ_opt\[ -4.176810958536834 \cdot 10^{-19} - 2.7755575615628914 \cdot 10^{-17}x_{2} - 3.122502256758253 \cdot 10^{-17}x_{1} + 0.0015948166443607303x_{2}^{2} - 0.001748548186795075x_{1}x_{2} - 0.0005642393727047848x_{1}^{2} - 0.004068980742835293x_{2}^{3} - 0.00889705628280991x_{1}x_{2}^{2} + 0.031008121446712326x_{1}^{2}x_{2} - 0.00462609155683142x_{1}^{3} - 0.13854894702932397x_{2}^{4} + 0.23890192550821426x_{1}x_{2}^{3} - 0.20074288167828302x_{1}^{2}x_{2}^{2} + 0.15886473368982973x_{1}^{3}x_{2} + 0.052300384255258056x_{1}^{4} + 0.21015195286319743x_{2}^{5} + 0.23612833587408735x_{1}x_{2}^{4} - 1.771046763880159x_{1}^{2}x_{2}^{3} + 2.325418396012682x_{1}^{3}x_{2}^{2} - 2.14073242822923x_{1}^{4}x_{2} + 0.13263464829810995x_{1}^{5} + 3.7316697440905724x_{2}^{6} - 8.975646042650421x_{1}x_{2}^{5} + 12.647030312088823x_{1}^{2}x_{2}^{4} - 13.157874488714592x_{1}^{3}x_{2}^{3} + 6.408472753225714x_{1}^{4}x_{2}^{2} - 4.70607484123265x_{1}^{5}x_{2} - 1.5648576613673593x_{1}^{6} - 3.1714641722949395x_{2}^{7} - 0.8235555675876345x_{1}x_{2}^{6} + 28.387225515889497x_{1}^{2}x_{2}^{5} - 77.62607234842116x_{1}^{3}x_{2}^{4} + 116.37607782247242x_{1}^{4}x_{2}^{3} - 80.11035758835642x_{1}^{5}x_{2}^{2} + 44.47069285639724x_{1}^{6}x_{2} + 0.0726000151063359x_{1}^{7} - 39.8533682575852x_{2}^{8} + 125.31127500418219x_{1}x_{2}^{7} - 220.58401873093374x_{1}^{2}x_{2}^{6} + 293.7720621346866x_{1}^{3}x_{2}^{5} - 266.9900870830618x_{1}^{4}x_{2}^{4} + 209.8387489805005x_{1}^{5}x_{2}^{3} - 59.10991674302298x_{1}^{6}x_{2}^{2} + 52.436849442791065x_{1}^{7}x_{2} + 18.20748677043684x_{1}^{8} + 14.506630503210488x_{2}^{9} - 9.440411382923457x_{1}x_{2}^{8} - 160.47240428813274x_{1}^{2}x_{2}^{7} + 747.8595775239182x_{1}^{3}x_{2}^{6} - 1695.7112578096178x_{1}^{4}x_{2}^{5} + 2153.382772341556x_{1}^{5}x_{2}^{4} - 1819.9767034601814x_{1}^{6}x_{2}^{3} + 855.579875742001x_{1}^{7}x_{2}^{2} - 351.2840811455434x_{1}^{8}x_{2} - 24.416052724567457x_{1}^{9} + 147.774293218062x_{2}^{10} - 583.0053139713531x_{1}x_{2}^{9} + 1169.7962341650689x_{1}^{2}x_{2}^{8} - 1923.8367515603788x_{1}^{3}x_{2}^{7} + 2237.5803557678155x_{1}^{4}x_{2}^{6} - 1909.954749738621x_{1}^{5}x_{2}^{5} + 974.4632412102649x_{1}^{6}x_{2}^{4} - 629.2192630930858x_{1}^{7}x_{2}^{3} - 33.20120287464316x_{1}^{8}x_{2}^{2} - 147.44786948201514x_{1}^{9}x_{2} - 65.28867531870668x_{1}^{10} + 447.63909430216984x_{1}^{2}x_{2}^{9} - 1691.4500588977705x_{1}^{3}x_{2}^{8} + 6377.2399025187315x_{1}^{4}x_{2}^{7} - 13592.058329862253x_{1}^{5}x_{2}^{6} + 15982.087567872013x_{1}^{6}x_{2}^{5} - 12839.154601026754x_{1}^{7}x_{2}^{4} + 7717.973145614064x_{1}^{8}x_{2}^{3} - 2808.720811515523x_{1}^{9}x_{2}^{2} + 900.4635509762891x_{1}^{10}x_{2} + 154.41745497532926x_{1}^{11} - 232.77539083718503x_{1}^{3}x_{2}^{9} + 250.74784828597632x_{1}^{4}x_{2}^{8} - 1604.7710599708485x_{1}^{5}x_{2}^{7} + 5427.010771019974x_{1}^{6}x_{2}^{6} - 7584.898114887506x_{1}^{7}x_{2}^{5} + 6647.784025839779x_{1}^{8}x_{2}^{4} - 4266.449127496643x_{1}^{9}x_{2}^{3} + 1643.2270736360429x_{1}^{10}x_{2}^{2} - 554.7169315237667x_{1}^{11}x_{2} - 74.66276382537889x_{1}^{12} \]
We plot its level sets inside the disk where the residual is constrained. Values outside this disk are extrapolations and are excluded from both the plot and the error calculation below.
using Plots
x1s = x2s = range(-r, r, length = 40)
inside = [u^2 + v^2 ≤ r^2 for v in x2s, u in x1s]
ϕs = [
inside[i, j] ? ϕ_opt(x => [u, v]) : NaN for
(i, v) in enumerate(x2s), (j, u) in enumerate(x1s)
]
contour(
x1s,
x2s,
ϕs;
levels = [-0.2, -0.1, 0, 0.1, 0.2, 0.5, 0.7],
aspect_ratio = :equal,
xlabel = "x₁",
ylabel = "x₂",
)
θ = range(0, 2π, length = 100)
plot!(r * cos.(θ), r * sin.(θ), label = "")
scatter!([0], [0], label = "")
We can check the Koopman identity using the finite-time Laplace average of ϕ_opt. An exact eigenfunction satisfies exp(-λ * t) * ϕ(S(t, x)) = ϕ(x) along a trajectory. Solve each trajectory once, then integrate using its dense interpolant. Tight tolerances matter because exp(-λ * t) amplifies integration errors.
import OrdinaryDiffEqTsit5
using QuadGK
function laplace_average(f, λ, x1, x2, T = 10.0)
prob = OrdinaryDiffEqTsit5.ODEProblem(
(v, p, t) -> [g(x => v) for g in F],
[x1, x2],
(0.0, T),
)
traj = OrdinaryDiffEqTsit5.solve(
prob,
OrdinaryDiffEqTsit5.Tsit5();
reltol = 1e-11,
abstol = 1e-13,
)
v, _ = quadgk(t -> f(traj(t)) * exp(-λ * t), 0.0, T; rtol = 1e-7)
return v / T
end
lap(x1, x2) = laplace_average(v -> ϕ_opt(x => v), λ, x1, x2)
laplace = [
inside[i, j] ? lap(u, v) : NaN for
(i, v) in enumerate(x2s), (j, u) in enumerate(x1s)
]40×40 Matrix{Float64}:
NaN NaN NaN NaN NaN … NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN … NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
NaN NaN NaN NaN -0.110708 -0.0401727 NaN NaN NaN NaN
NaN NaN NaN NaN -0.111655 -0.0455329 NaN NaN NaN NaN
⋮ ⋱ ⋮
NaN NaN NaN NaN 0.210627 0.273319 NaN NaN NaN NaN
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN … NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaNWe report the maximum absolute discrepancy on the disk. Unlike a matrix norm over the enclosing square, this error excludes unconstrained points and does not grow merely because the grid is refined. It checks the eigenfunction identity over a finite time, rather than providing a bound on the error with respect to the exact eigenfunction.
maximum(abs, laplace[inside] - ϕs[inside])0.022892164571144027The disk is not forward invariant: some trajectories starting inside it leave it, so a small residual on the disk does not guarantee a small Laplace-average discrepancy for all those initial states. Increasing T does not fix extrapolation outside the approximation domain. For initial states in the smaller disk of radius r / 2, the discrepancy is much smaller:
inner = [u^2 + v^2 ≤ (r / 2)^2 for v in x2s, u in x1s]
maximum(abs, laplace[inner] - ϕs[inner])6.191118991377698e-5The finite-time identity discrepancies above need not match the error reported in the book, especially when eigenfunctions have different scalar normalizations.
This page was generated using Literate.jl.